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INTRODUCTION 


Advances in computational techniques and computers in the past few years make it 
practical to compute the steady inviscid flow field about complex three-dimensional 
bodies, such as the Space Shuttle orbiter or other advanced entry vehicles, in their 
actual supersonic or hypersonic flight environment. The inviscid flow field provides 
surface pressures, which can be integrated to obtain aerodynamic loads, and other 
flow properties which are required to calculate surface heating rates (ref, 1) needed 
to define the thermal environment. 

These vehicles enter the atmosphere at relatively large angles of attack which 
will lead to one of two classes of problems. (See fig. 1.) If the angle of attack 
is moderate (25® to 35® for Shuttle-like vehicles), the subsonic portion of the flow 
field is generally confined to the vehicle nose. Several papers (refs. 2 to 5) have 
presented time-asymptotic methods of efficiently solving the three-dimensional invis- 
cid flow over blunt-nosed bodies at moderate angle of attack where the subsonic 
region is relatively small. These solutions provide a data surface, downstream of 
the subsonic region, on which the local flow velocity is supersonic. Several papers 
(refs. 6 to 9) have presented methods for continuing the solution downstream in the 
supersonic region by using spatial marching techniques. Since these techniques use 
spatial marching, they require relatively low computer storage. These methods have 
been shown to provide good results unless additional embedded pockets of subsonic 
flow are encountered (such as near the leading edge of wings) where the spatial 
marching techniques break down. 

When the angle of attack is large (greater than 35°), the subsonic region is no 
longer confined to the nose but extends much farther downstream and can envelop much 
of the lower surface. (See fig. 1.) Since the entire subsonic region must be com- 
puted simultaneously, the time-asymptotic portion of the solution will require many 
more grid points than required for the moderate angle-of-attack case discussed pre- 
viously. Therefore, codes must be structured for computers that have the storage and 
computational speed required to solve this type of problem (which may require 90 000 
or more grid points). Thus, existing time-asymptotic methods (refs. 2 to 5), either 
because they are constructed for scalar computation or because of the coordinate 
system used in the code, are not suited to solve the flow over complete vehicles with 
complex three-dimensional geometries at high angle of attack. In reference 10, it is 
shown that a vector-processing computer is ideally suited for solving this type of 
large flow-field problem. 

The present paper presents a time-asymptotic method that is being developed for 
the CDC® CYBER 203 vector-processing computer which will be able to solve the flow 
over complex three-dimensional bodies (Shuttle-like geometries at large angle of 
attack) where large embedded subsonic regions occur. Results are presented in this 
paper which demonstrate the capability of the code HALIS (J^igh Alp ha ^nviscid 
_^olution) to compute the flow field over the Space Shuttle orbiter at large angles of 
attack. Additional applications of the HALIS code may be found in reference 11 for 
comparisons with Shuttle tunnel data and reference 12 for comparison with Shuttle 
flight data. The components of the Euler equations as well as the development of the 
transformed equations in both the spherical and cylindrical coordinate systems are 
presented in appendix A. The Euler equations in the spherical coordinate system 
along the ray 0 = n are given in appendix B. A stability analysis is presented in 


appendix C by M, J. Hamilton. The method for determining the shock velocity and 
other flow properties at the shock wave is outlined in appendix D. 



SYMBOLS 

local speed of sound, 
pressure coefficient 
total energy, 
internal energy, e/Uoo 

dummy variable used in generalized description of derivative forms 
total enthalpy, H/U^^ 
static e n tha Ipy , h / 

total length of vehicle, 1290 inches to hinge line 
reference length, ft 
Mach number 

outer unit normal to body 
outer unit normal to shock 
pressure, P/p„U^2 
gas constant, f t^/sec^-®R 
cylindrical coordinate, in. 
spherical coordinate, in. 
entropy, S/R 

vector defined by equation (11) 
time, t/(L^/U^) 
velocity, ft/sec 
r-velocity component, 
f -velocity component, 
z -velocity component, u^/U^ 

0-velocity component, Ug/U^ 

((rve loci ty component , u 
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^velocity component, 

V total velocity vector, V/U^ 

A# 

Vg shock velocity, V^/U^ 

W,F,G,H,Q vectors defined in equations (A8) 

W,P,G,H,Q vectors defined in equations (A2) 

W,F,G,S,Q vectors defined in equations (B4) 

z cylindrical coordinate 

a angle of attack, deg 

r distance between body and shock, defined in figure 7 

Y ratio of specific heats 

e smoothing coefficient 

r\ transformed cylindrical coordinate, equations (A9) 

rj transformed spherical coordinate, equations (A3) 

9 spherical coordinate 

^ transformed cylindrical coordinate, equations (A9) 

p density, p/p„ 

(|> cylindrical coordinate 

(j> spherical coordinate 

(i» transformed cylindrical coordinate, equations (a9) 

(\) transformed spherical coordinate, equations (A3) 

0 ) transformed spherical coordinate, equations (A3) 

Subscripts: 

b body surface 

max maximum 

s shock surface 

w wall surface 

00 free-stream conditions 
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Notation over symbols: 

~ dimensional variable 

vector 

^ unit vector 

spherical coordinate system 
= 4) == 71 plane in spherical coordinate system 


CCXDRDINATE SYSTEMS 

The method of solution presented in this report requires that the outer boundary 
of the solution space correspond to the bow shock which envelops the entire vehicle 
and that the inner boundary correspond to the vehicle surface. Also, it is required 
that all coordinate lines extending between the body surface and the shock intersect 
the body surface only one time. When the vehicle to be represented is a short blunt 
body, a spherical coordinate system naturally satisfies these requirements. However, 
when the body is blunt and many nose radii in length, the use of a spherical coordi- 
nate system would result in a highly skewed physical grid. In order to avoid this 
problem and still satisfy these grid requirements, the physical grid is constructed 
by combining a spherical with a cylindrical coordinate system. Such a combined 
coordinate system is shown in figure 2. In this right-hand system, the spherical and 
cylindrical coordinate system are coupled at the plane 0 = n/2, z = 0 where the 
two systems are coincident. 

In the physical domain, computations are made in the spherical system over 


^<♦<1 

■? < 0 < n 
^ J 

and in the cylindrical system over 


r, < r < r 
b s 


^<♦<1 ) 


0 < z < z 


max 




(2) 


A typical representation of this physical grid system is shown in figure 3(a), a 
symmetry-plane view, and in figure 3(b), a cross-flow-plane view which corresponds to 
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the exit plane of figure 3(a), in order to improve the clarity of these two figures, 
a number of grid lines have been removed. As can be seen, the physical coordinate 
system is skewed and nonorthogonal. 


GEOMETRY DESCRIPTION 

The vehicle configurations used for making flow-field calculations with the 
present code are first modeled with QUICK - a geometry program described in refer- 
ence 13. This model provides a smooth analytic description of the vehicle geometry 
in a "local polar coordinate system," r^ = r^(z,(|)), with continuous surface deriva- 
tives over the vehicle, Hie appropriate subroutines from QUICK have been included in 
the present flow-field code so that the geometry models generated by QUICK can be 
used in the present code for making flow-field calculations. 

Two versions of the geometry of the Space Shuttle orbiter (shown in fig, 4) have 
been modeled with QUICK as shown in figure 5* The first (fig. 5(a)) is a reasonably 
accurate model of the actual vehicle with the canopy, vertical tail, and orbital 
maneuverable system (QMS) pod removed. The second (fig, 5(b)) has the same lower 
surface shape and the same profile for the upper symmetry plane as the first, but the 
region between the leading edge of the strake or wing and the upper symmetry plane 
have been replaced with elliptical segments as shown by dashed lines in section A-A 
of figure 4. This process simplifies the leeside geometry and makes calculations in 
this region easier but does not affect the results obtained on the windward side 
since the cross-flow velocity goes supersonic near the wing tip. Since the flow in 
the lower part of the shock layer on the leeside of a vehicle at large angle of 
attack is viscous dominated (ref, 14) and cannot be accurately modeled with the 
inviscid equations of motion, alteration of the inviscid solution on the leeside (due 
to the modification of the leeside geometry) is of little consequence. Thus, all the 
results for the flow field over the Space Shuttle orbiter presented in the present 
paper were obtained with the modified geometry shown in figure 5(b), 


METHOD OF SOLUTION 
Flow-Field Equations 

The flow field of interest in this paper, will contain one or more internal 
shock waves. Specifically, a cross-flow shock will be located on the leeside of the 
vehicle near the upper symmetry plane, and under certain conditions a strake/wing 
shock may appear. These shock waves must be properly resolved; thus, we have chosen 
to "capture" these shocks since this technique is compatible with the vector- 
processing characteristics of the CYBER 203 computer. Thus, in this paper, the time- 
dependent, three-dimensional, compressible Euler equations are integrated in the weak 
conservation form. These equations can be written in the general form in the spher- 
ical coordinate system as 


aw 

at 


+ 


a? 




(3) 


where the vectors W, F, G, and H represent the usual conserved quantities in the 
spherical coordinate system and the vector Q contains all the terms that arise from 
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the use of a non-Cartesian coordinate system. In appendix A, the form of these vec- 
tors is given for spherical coordinates in equations (A2), Likewise, the equations 
can be written in the cylindrical coordinate system as 


^ ^ 

at ar a 4) 



(4) 


where the vectors W, F, H, G, and Q are given in equations (A8) . 

The flow equations (eqs, (3) and (4)) are transformed from the physical domain 
to separate computational cubes, one corresponding to each coordinate system, by the 
following equations which allow the description of a general body in terms of its 
radius. For the spherical system. 


r - 

h = *1 1 1 — :: — 

r^(t,4),0) - 


(0 < n < 1) 


^ = UuiZi 

It 


(0 < ^ < 1 ) > 


(5) 


w 


Tt - 0 
lt/2 


(0 < 0) < 1) J 


and for the cylindrical system. 


r] = 


4> = 


5 = 


r - r ( ()),z) 
b 


r^(t, <|),z) - r^^{ 4),z) 


<1> + u/2 
7t 


z 


z 

max 


(0 < T) < 1) 

(0 < 4 > < 1 ) \ 


(0 < 5 < 1) J 


(6) 


Details of the transformation are outlined in appendix A, Although the physical 
domain is transformed to a computational domain, the velocity components are not 
transformed, and they retain the same magnitude and direction in the computational 
mesh as they did in the physical grid. In this paper, all computations are made by 
using equally spaced meshes. 
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Numerical Procedure 


The computational data base is arranged as shown in figure 6 . The data base is 
arranged in an interleaved manner (ref. 15). Interleaving, the sequential storage of 
the data base associated with each computational plane, allows the efficient use of 
the virtual memory system of the computer should the size of the code ever exceed the 
available central memory. 

Interior grid points .- The following discussion applies to all interior points 
in the computational space except for those lying in the plane 9 = u, which will be 
treated in a separate section. The governing equations (eqs. (A 6 ) (spherical) and 
(A12) (cylindrical)) are integrated in time with an unsplit MacCormack differencing 
scheme (ref. 16). In the predictor step, the plane-by-plane integration is carried 
out from the plane 9 = u to the plane z = with forward differences in each 

direction. (See fig. 6 .) 

The plane where the two coordinate systems are coincident is treated as part of 
the spherical system when integrating the governing equations. To construct the for- 
ward streamwise derivatives in the predictor step at the plane requires that a plane 
of data be established at 0 = 3. ^ 9 , This temporary plane of information is 
obtained by interpolation from the planes of data in the cylindrical system immedi- 
ately downstream of the plane z = 0. The thermodynamic properties transfer directly 
from one coordinate system to the next; however, the interpolated velocity components 
in the cylindrical system must be mapped back to the appropriate velocity components 
in the spherical system. 

To preserve the interleaving concept, the corrector step is begun in the plane 
z == 2j^ax sweeps back to the plane 9 = by using backward differences. Since 

backward differences are used, no interpolation is needed when computing at the plane 
9 = n/2 in the spherical system. Also, backward derivatives for use in the first 
computed plane in cylindrical system can be constructed directly from information in 
the plane 9 = n/2 since it is part of both coordinate systems. This construction 
must take into account that in the plane 6 =* ii/2, 


u 

z 


-u 
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Differencing at boundaries .- The governing equations must be integrated at both 
the vehicle surface and the bow shock. To preserve the second-order spatial accuracy 
of the T) derivatives at these boundaries, the method of Abbett (ref. 17) is fol- 
lowed and these derivatives in the spherical system are defined in the predictor step 
as 


5f 



af 


5t1 


s 


— (f - f ) 

- w +1 w 

Ati 


S S-1 

Ati 




(7a) 
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and in the corrector step as 


5ti 


w 


= J_(-2f + 3f , - f 

- w w+1 w+2 

At) 


— =— (2f - 3f + f J 

- *- s s-1 S.-2 

5ti ^ Ati 




(7b) 


Similar equations are used in the cylindrical system by replacing t) with t). As 
shown by Abbett, the net effect of using these derivatives in the MacCormack differ- 
encing scheme is to have applied a second-order accurate, three-point backward dif- 
ference at the boundary. 

Shock derivatives.- Derivatives on the shock surface are necessary to update the 
physical grid and to determine post-shock properties. When evaluating the streamwise 
shock derivatives 5r /90 and 5r /bz, it is necessary to use noncentered four-point 

s' s' 

derivatives of the form 


f . 

1 


1 

6 A 


^^2 


6f. ^ 
1-1 


+ 3f . + 2f . ) 

1 1+1 


to prevent downstream disturbances from moving upstream along the shock and causing 
oscillations in the shock about the stagnation point. At the outflow boundary a 
three-point backward difference derivative is used for Br^/Sz . in the plane 0 = 
which at the shock is really a point and not a surface, best results are obtained 
when a centered two-point derivative is used for br^/QG. Formulation of the stream- 
wise shock derivatives by using noncentered four-point derivatives about the juncture 
of the two coordinate systems requires that interpolated values of the shock position 
be determined in both coordinate systems for both the predictor and corrector steps. 

in the cross-flow planes, dr /dtf) and dr /d<l> are approximated by two-point 
centered derivatives. 

Plane 0 = The governing equations in the spherical coordinate system 
(eq. (A1 ) ) are singular at 0 = -ji. in appendix B, a reduced set of equations, valid 
at 0 = 7 t, are developed. The desire is to integrate these equations in a manner 
that is consistent with the integration scheme, that is, unsplit MacCormack scheme. 
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used in the rest of the symmetry plane. (See appendix B.) The equation to be 
integrated is equation (B6), which is restated below in the computational 
coordinates; 


at 




+ + + + 5 . ^ . Ad = n 

at ae ae ao) ^ ^ aft ?w.i - - ) ^ afi -I - - ) ^ 


at) ar 5 t| 


ari 


ae 


ao)^ 


a4>> 


ae 


( 8 ) 


ar)\a(t> b(\>/ 


The integration of equation (8) requires information from the first two data 
planes shown in figure 6. These two planes are shown in greater detail in the fol- 
lowing sketch: 




In appendix B, the number of equations represented by equation (8) was reduced 
from five to four (at 0 = it) by restricting the integration of equation (8) to 
points in the symmetry plane. To further simplify the integration, computations are 
done in the windward symmetry plane, ^ = -%/2» In equation (8), the integration 
with respect to t is identical to that used in the rest of the flow field. In the 
predictor step, the derivative 5G/5a> is taken to be 


1 

Ao) 


5(, 


9=7C-A9f (j)=-“f 




(9a) 


and in the corrector step, the derivative 55/do) is approximated by 




(9b) 


When forming the quantities G ^0=7t-A0f the change in sign of Uq must be 

taken into account when going from $ = -it/2 to $ = n/2. (See eq. (B2a).) 


9 


The formation of the second derivatives in equation ( 8 ) can be simplified by 
examining the behavior of ^derivatives with respect to ^ in the physical plane* 
First take 5S/5^ with S = fu- so that 


a(fu-) au- 

^ = u ^ + f 

9(f) ^ 9$ 94) 


In the definition of U 7 on 6 = u (eq* (B2b)), 

4> 


i = -U 




^-n/2 


sin ^ 


or 


&u- 


d(t> 


= u 


(t>=-u/2 


(f>=-n/2 


along the line ^ = -n/2. Finally, 




(})=-ll/2 


f + u- 


M 

a^ 


Now S (eq. (B4d) ) can be written as 


S = u- T 


where 


rp 

rpu- 

rpuQ 
r(E + p) 


(10) 


( 11 ) 
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The substitution of T for the dummy function f in equation (10), along with the 
fact that T is an even function about the symmetry plane, leads to the expression 


5^ 


<|f=-7c/2 


T 


From equations (A4), 


5 _ ^ ^ ^ 

5^ 9$ 5t) 5^ 


and it has been established that 


^ = 
9^ 


0 


on 


0 = 7 t; thus. 


OS _ 94^ 9S _ ^ 

5? dt ait ^ 




This allows us to write the term 


ari\a(t) aity 


found in equation (8) as 


^ a 


<Jf=-7i/2 


Thus, the differentiation can be handled the same as for the other rj derivatives. 
In a similar manner, in the predictor step. 


au) a /a<t as '' 


Aw a 0 




_ 

(?Ug) 


\ai 

ait/ 

&=it-A0 

0=lt 

(J)=-lt/2 


(12a) 


1 1 


and in the corrector step 


^ A. M ^ 

a?y 



— 

/ _ 5u_\ 

* 

Ao) 69 


. p f 1 a \ 


e=ir \a^ d<t> / 

0=ii-A0 






(12b) 


Again, care must be taken to properly account for the change in sign of u^ in the 
symmetry plane ^ = it/2. 

T?he wall boundary conditions, in the plane Q = % are handled as in the rest of 
the flow field* Procedures for computations on the shock boundary are the same as 
those used in the rest of the field although some terms in the equations must be 
altered because of the coordinate singularity. These changes are included in the 
discussion of the shock computations in appendix D* 


Boundary Conditions 

Wall boundary *- For inviscid flow, the wall boundary condition requires that the 
velocity component normal to the wall be zero and that the surface entropy be con- 
stant at the post-normal-shock condition. In addition, in the steady state, the 
total enthalpy at every point in the flow field must be a constant. 

The following procedure is used to satisfy these boundary conditions. The con- 
tinuity and momentum equations are integrated along the surface. The computed wall 
density and the constant wall entropy condition are used to determine the internal 
energy and the pressure at the surface. Since, in general, the physical grid is not 
normal to the body surface nor the cross-flow and axial components of the velocity 
tangent to the body surface, the surface velocity boundary condition becomes 
^b * ^b ~ where in the spherical system. 


V, = u- , e- + u- , e- + u. ^e. 
b r,b r ^ 0,b 0 


(13) 


is the total velocity vector at the surface and n^ is the outer unit normal at the 
wall. In general, the velocity vector formed from the velocities determined by 

integrating .the momentum equations at the wall does not satisfy the tangency condi- 
tion = 0 and lies outside the surface tangent plane. To correct the total 

velocity vector at the wall so that it lies in the surface tangency plane, a scheme 
similar to that outlined by Kutler (ref. 18) for use in a three-dimensional marching 
code is used. Unlike Kutler, the computed thermodynamic properties at the wall are 
assumed to be correct, and a set of surface velocities are sought which are con- 
sistent with those thermodynamic properties and the specified wall velocity condi- 
tion. Although the resulting wall properties are inconsistent at any instant in 
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time, the solution at the wall converges as the global solution converges. For the 
general case. 


b,N 




( 14 ) 


where 


= [2(H - (15) 

b w 

is the total velocity at the surface, Vj^ ^ is the corrected velocity vector on the 
body surface which lies in the surface tangent plane and 


n^ = 


► A A 

w - (V. • n. )n, 

b,o b,o D t 


V. 


A A 

^ - (V • a )n 
),0 b,o D I 


(16) 


is the unit vector in the surface tangency plane where V, 
ity vector at the wall. 


is the computed veloc- 


In the spherical system, the outer unit normal at the wall is given by 


% = 



5f 


b ^ 
— e- 


r sin 9 5^ 


1 


^r sin 9 5(|) 



! 1/2 


(17) 


By combining equations (13) to (17), the following expressions can be determined for 
the velocity components at the surface: 


^^r,N ^byr,o ^ 2/1 


(18a) 


Di dr^^. 

“e,N " ^b|“e,o ’^1/^3 

^b 2 


(18b) 


°1 \^^b 


(18c) 
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where 


c , U- . ^ ^ ^ 

^ r r sin 0 


ftr ^^2 

, 1 

°2 = ' ^ 


“1 2 


6r, 


3(1) 


^ (r sin 0) 

D 


D \2 


ar. 


n 2 


^^0,0 °i -wr^h^2'> 


3r, 


—1 9 ^ V 2 


"4>'° ■" ■;7r"b°2 

0(p 


In the cylincJrical system, a similar analysis leatis to the following expressions for 
the surface velocity components: 


'^r,N ‘^b|^r,o ” ^ 2)1 ^ 


(19a) 


'^4),N = ^b 


arb/ 

"♦,0 ■" °i 8r/‘V2> 


(19b) 


I 

^,N =%hz,o °1 -^I^2p3 


(19c) 


where 




ar. 


= U ~ U 


r, - u 


1 r,o (|),o\6(t) lb/ z 9z 


/ ar \2 ar 

D = 1 + —±r + — ^ 

2 \a(() hi dz 
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2'^V2 


D 


3 



r,o 



< t >/0 1 


9(t) 


“1 2 





+ D. 



Shock boundary The bow shock, which is the outer boundary of the flow field, 
requires a time-dependent boundary condition, since the post-shock properties and 
shock location are a function of the computational time increment. The details of 
how the bow-shock boundary is handled in both coordinate systems are given in 
appendix C, Briefly, the post-shock properties are determined by the post-shock 
pressure and the local inclination of the shock to the free-stream velocity vector. 
The post-shock pressure is taken to be the pressure at the shock location determined 
by the integration of the flow-field equations. From this pressure and the shock 
geometry, new post-shock conditions and a shock velocity can be determined. If 
incremental movements of the shock are assumed to be small, then the change in radial 
location of the shock can be written in the spherical system as 


= Vg dt 


(20a) 


and in the cylindrical system as 


= Vg dt 


(20b) 


A MacCormack scheme is used to integrate these equations to update the position of 
the shock. At convergence, the computed post-shock pressure will give a zero shock 
speed and thus a stationary shock wave. 

Outflow boundary ,- Computations are always carried to a point where the axial 
flow at all grid points in the cross-flow plane is supersonic. We then integrate the 
regular interior point equations in the outflow plane by approximating axial deriva- 
tives with two-point backward differences in both the predictor and corrector steps. 

Smoothing function ,- Historically, flow-field calculations, with the MacCormack 
integration scheme, have required the addition of some damping to maintain numerical 
stability; the HALIS code is no exception. We have chosen to use the fourth-order 
smoother proposed by Barnwell (ref, 19). The smoothing is done on data in the compu- 
tational plane and is applied at each time step after completion of the corrector 
step. Because of CPU (central processing unit) storage considerations, the primitive 
rather than conserved variables are smoothed; this has proven to be satisfactory. 

The smoothing is formulated as follows in the spherical system: 


^t+At ^ ^t+At 
CO Tl/(^/C0 


00 ) 


5^ 


,.->4 

(An) 

an j 


( 21 ) 
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The same equation holds in the cylindrical system when ^ is substituted for co, 
for and t) for r\* 

Ihe fourth-order derivatives are formulated by using regular five-point central 
differences. In the n and rj directions^ at one point off the wall or shock, 
third-order derivatives using four-point differences are substituted for the fourth- 
order derivatives in equation (21). No damping is used in the n or t] directions 
at the wall or the shock. A third-order derivative is substituted for the fourth- 
order 5 derivative in the next to the last plane in the cylindrical system whereas 
no smoothing is applied in the ^ direction on the outflow plane. Barnwell 
(ref. 19) has shown that 0 < e < 1/24. A value of e of 0.025 has been found to 
provide adequate damping without distorting the computed flow field. 

Stability analysis .- A closed-form solution for the MacCormack differencing 
scheme when applied to the three-dimensional Euler equations is very difficult, if 
not impossible, to obtain. In appendix C, a rigorous derivation of the stability 
criterion for the Brailovskaya differencing scheme is given. It is further shown 
that under "worst-case" assumptions, stability for the MacCormack scheme reduces to 
that of Brailovskaya. Thus, the stability criterion used in this work is based on 
equation (C24) • 

At each point in the computational grid, the time integration is allowed, to 
advance at its own local allowable time step. However, the streamwise distribution 
of time step should be smooth at the juncture of the two coordinate systems. To 
guarantee this smoothness, the magnitude of the local time steps in the spherical 
system is adjusted within the limits of the stability criterion. 

The following usual definition of the limiting time step is used: 


At 


1 



( 22 ) 


where maximum eigenvalue occurring in the stability matrix. In the 

spherical coordinate system. 



(23) 
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where a is the local speed of sound and 



(rg - rj^) An 


5r 5r 

(1 - n) — ^ + n — 


_J 

n Acp 


In the cylindrical coordinate system. 
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h 


1 


(1 




n A(|; 


(r 

s 





ar 


f = 


(r^ - An 


(1 - Tl) 


9z 


+ T) 


9r 

c 

9z 


z 

max 


Further analysis indicates that the value of the smoothing coefficient used will 
affect the magnitude of the allowable time step size. From this analysis, equa- 
tion (22) is replaced with the following equation when determining the maximum 
allowable local time step: 


At < 




96e + (1 + 192e) 




1/2 




max 


Thermodynamics For perfect gas computations, the pressure is computed by 
using the state equation p = (y - the speed of sound from the equation 

= Y(Y - 1)e. A modified HALIS code has been written which incorporates 
equilibrixim-air chemistry. Hiis version uses the equilibrium-air curve fits of 
P = p(pf®) ^ = a(p,e) done by Tannehill and Mugge (ref. 20). Some difficulty 

with convergence can be experienced when using these curve fits because of slight 
differences in thermodynamic properties at the juncture of fitted-curve segments. 

This problem has been overcome by checking for discontinuities in the thermodynamic 
properties at the juncture of fitted-curve segments and, where they exist, generating 
a smooth transition from one curve segment to the next. To reduce computation time 
when running this code, the nonvectorized chemistry code is called every 25th itera- 
tion. In the interim, a local effective y used in computing the pressure. This 

procedure keeps the time for equilibrium chemistry computations to approximately 
20 percent greater than for a perfect gas. 


Initialization 

This method of solution, which is posed as an initial value problem, has 
required the development of an initialization procedure for the flow field about a 
complex three-dimensional geometry. In the following text, the flow initialization 
procedure as currently used in the HALIS code is described and is the result of the 
investigation of numerous initialization procedures for different parts of the flow 
field. 


Surface properties .- Starting from the stagnation point, the Newtonian pressure 
on the surface is determined in the symmetry plane. On the leeside, the pressure is 
only allowed to drop to six times the free-stream pressure. The multiple six has no 
physical significance, but, using smaller values has, in some cases, led to negative 
pressures on the leeside. Then, the meridional pressure distribution is set in each 
plane by a sine function. Once the pressure distribution is set, an isentropic 
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expansion from the stagnation point gives the density and internal energy on the 
surface. Ttie total velocity on the surface is then computed by assuming that the 
total enthalpy remains constant throughout the flow. To split the total velocity 
vector into its three components, the simplified streamline method of DeJarnette 
(ref. 21) is used to determine the direction cosines of the streamline at a point on 
the surface. Each component of velocity is then determined from the local total 
velocity and the appropriate direction cosine. In addition, this procedure auto- 
matically satisfies the surface boundary condition v • n =0. 

b 

Shock configuration and properties ,- As described previously, the bow shock 
serves as the outer boundary of the coordinate system. We have found that it is 
imperative that a reasonable guess be made of the initial shock shape, that is, that 
the shock not be required to make large spatial adjustments. Initially, rotated 
paraboloids were used as initial shock shapes; this works well for short bodies at 
moderate angles of attack. However, this procedure breaks down for long bodies at 
high angles of attack as the shock on the windward side of the body will intersect 
the body surface. 

Currently, a segmented approach is used to determine an initial shock shape. In 
reference to figure 7, the shock shape in the symmetry plane is broken into four 
separate curve segments, A-B, B-C, C-D, and D=E.. Curve B-C is the locus of points 
which are equidistant from the body surface when measured along a body normal. The 
separation distance is the shock standoff distance, determined from the relation 
(ref. 22) 


6 

s 


0.78 - 

P /P ^ 


Curve C-D is a parabola which matches the slope and position of curve B-C at 0 = jt. 
The point of intersection of the curve C-D and the ray 0 = n/2, $ = n/2 is a vari- 

able r, which depends on the angle of attack. Curve D-E is a straight line that 
matches the slope and position of curve C-D at F; thus, the choice of F determines 
the shape of the leeside shock. Curve A-B is a second-order curve which matches the 
slope and position of curve B-C at 0 = 'jz/2, $ = -ti/ 2 and has a slope equal to the 
body slope at z = In cross section, the shock shape is taken to be a circle 

whose diameter is the distance between the shock in lower and upper symmetry planes. 
Typical initial shock-shape cross sections are shown in figure 8 for several loca- 
tions of z. At first glance, this procedure appears very cumbersome. However, it 
is easily automated and is only dependent on picking a reasonable value of T, which 
for the specific geometry used in this paper varied from 85 inches to 140 inches over 
angles of attack from 25° to 45°. Once the shock geometry is established, the same 
procedure as outlined in appendix D for handling the unsteady shock wave can be used 
to find the initial post-shock properties by using steady-state shock relationships 
obtained by setting the shock velocity to zero in the relations of equations (D21), 
(D22), and (D23) . 

Interior grid points .- Ihe interior grid points are initialized by first forming 
the conserved quantities W, W, and W, as defined in appendixes A and B, along the 
shock and wall and then determining the values at the interior nodes by linear inter- 
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polation, Ihe primitive variables are then separated from the conserved quantities 
and the pressure at each node is defined from the thermodynamic relation 


P = P(p»e) 


( 26 ) 


HALTS configuration and operational statistics .- The requirements on grid 
resolution and the CYBER 203 CPU memory size limit the extent of the modified Space 
Shuttle geometry (fig. 5(b)), which can be modeled in the present HALTS code. These 
restraints limit the modeling to the first 650 inches, which represents approximately 
50 percent of the vehicle length measured from the nose to the body flap hinge 
line. Tn figure 6, each plane has 11 points in the r\ direction (distance between 
body and shock) and 32 points in the (\> (meridional) direction. There are a total 
of 90 planes in the streamwise direction, the first 15 in the co (or 0) direction 
in the spherical system and the rest in the ^ (or z) direction in the cylindrical 
system. Thus, the computational grid contains 31 680 grid points. The code requires 
that 22 variables be stored at each grid point. These variables include the old and 
new values of the flow-field variables, the physical coordinates, the local time 
step, and the smoothing variables. These storage requirements along with other stor- 
age needs and the associated computer code require all the CYBER 203 central memory 
of 10^ 64-bit words. Tn its present form, the code will generate converged solutions 
in 700 to 1500 global iterations for angles of attack from 20® to 45®. This corre- 
sponds to CPU time ranging from 2200 to 4500 seconds. A solution is considered to be 
converged when every point in the flow field meets the following criterion: 


t+At t 
P ■ P 
t 
P 


< 10 


-4 


(27) 


The HALTS code has a built-in restart capability which allows unconverged solutions 
to be continued or new solutions to be started from previously converged solutions. 
Using an old, converged solution to start a new one reduces by 10 to 15 percent the 
number of iterations required for convergence. 


RESULTS AND DISCUSSION 
Surface Boundary Conditions 

The surface boundary condition plays an important role in the development of the 
subsonic axial flow along the windward surface. Some authors (refs. 3 and 9) have 
avoided the subsonic flow problem by using solution techniques which do not directly 
impose one of the surface boundary conditions, namely that of constant surface 
entropy. Allowing the surface entropy to take on a value determined by integration 
of the flow-field equations, in general, leads to a surface entropy lower than that 
required by the constant entropy condition. This effect is illustrated in fig- 
ure 9(a) where center-line surface entropy is plotted as a function of Shuttle 
vehicle length for a constant and variable surface entropy condition. The lower 
surface entropy values in turn produce higher values of axial Mach number (fig. 9(b) ) 
but leave the surface pressure distribution unaltered (fig. 9(c)). Thus, it would 
appear that, if surface pressures are the only required product, computations may be 


20 


made by using a variable surface entropy. However, not imposing the constant wall 
entropy boundary condition will alter the distribution of flow variables through the 
shock layer. These altered distributions make such computed flow fields incompatible 
with three-dimensional and quasi -three -dimensional boundary-layer codes which rely on 
the solutions from inviscid flow-field codes for boundary-layer edge conditions. The 
extent to which these distributions are altered is shown in figure 10 where variable 
distributions through the shock layer at z/l = 0.1, 0.25, and 0.45 on the windward 
symmetry plane center line are shown for a HALIS computation of a flow field about 
the Shuttle vehicle at = 10.3 and a = 25®. Figure 10(a) is an entropy plot for 

both the case of a constant and variable wall entropy. Differences between the two 
solutions increase with increasing z/L, with these differences confined to the lower 
third of the shock layer. In fact, at z/L = 0.1 the solutions are virtually the 
same. Figures 10(b) and 10(c) show plots of u^r the axial component of velocity, 
and the density. The results are similar to those seen in the entropy plot. Fig- 
ure 10(d) is a pressure plot which shows that the different entropy boundary 
conditions have no effect on the pressure distribution through the shock layer at all 
three axial locations. In figure 11, plots similar to those shown in figure 10 are 
presented for HALIS computations of flow over the Shuttle vehicle at = 10.3 and 
a = 45®. Figures 10 and 11 show that differences between constant wall entropy and 
variable wall entropy are more pronounced at the higher angle of attack. These dif- 
ferences between the two solutions extend to the lower one-half of the shock layer 
at z/L = 0.45 for all three variables: S (fig. 11(a)), u^ (fig. 11(b)), and p 
(fig. 11(c)). However, in figure 11(d) the pressure distribution through the shock 
layer is unaltered by the choice of wall entropy boundary condition. 


Surface Axial Mach Number 

When the proper surface boundary conditions are applied, the area of the wind- 
ward surface having a subsonic axial velocity component should increase with angle of 
attack until the axial velocity component over the entire windward surface is sub- 
sonic as illustrated in figure 1 for the high-angle-of-attack case. Since the HALIS 
code requires a supersonic outflow boundary, it was necessary to alter the modified 
Shuttle geometry to include a slight expansion region aft of the 650-inch plane to 
raise the surface axial Mach number to a supersonic value. 

In figure 12, the windward-surface center-line axial Mach number distributions 
are shown for angles of attack between 25® and 45® and a free-stream Mach number of 
10.3. According to these results, if an initial data plane were established at 
z = 150 inches, then a spatial marching code such as STEIN (ref. 7) should be able to 
compute the flow over the vehicle at angles of attack up to 40®. However, practical 
experience with the STEIN code has shown that it is impractical and/or impossible to 
march into flows where M^ approaches 1. To determine, for this body and free- 
stream Mach number, what the limiting angle of attack for a STEIN solution would be, 
the HALIS code was used to establish an initial data plane at z = 1 50 inches. 
Complete STEIN solutions over the first 650 inches of the vehicle were obtained for 
angles of attack of 25®, 30®, and 35®. Center-line axial Mach numbers from STEIN 
solutions at a » 25®, 30®, and 35® are plotted at 100-inch increments in figure 12 
and agree with the complete HALIS solutions at the same angles of attack. For 
a = 40®, the STEIN solution failed after the first step away from the initial 
supersonic data plane. A broader view of how the subsonic axial Mach number region 
spreadr over the windward surface with increasing angle of attack can be found by 
comparing the parts of figure 13, which is a three-dimensional representation of the 
vehicle geometry being used with the physical grid superimposed on a black body. The 
white part of the vehicle surface represents the area of subsonic surface axial flow 
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and the black part represents the supersonic surface axial flow. The expansion 
region appended to the geometry can be seen at the aft end of the vehicle. At an 
angle of attack of 40°, a region of subsonic axial flow appears along the strake of 
the unmodified vehicle (fig, 5(a)), which appears as a chine along the modified 
vehicle. This subsonic axial flow is the result of the large expansion about this 
region at the higher angles of attack and the turning of the flow toward the upper 
symmetry plane which accompanies the expansion. At a = 42,5°, there are three 
regions of subsonic flow: the region from the stagnation point aft, the region about 

the strake, and the region near the aft end of the vehicle. This last region near 
the aft end of the vehicle is the result of the flow recompressing after it has 
expanded down the vehicle from the stagnation point. At an angle of attack of 45°, 
the axial flow over the entire windward surface is subsonic. 


Code Validation 

To validate the HALIS code for flow about complex three-dimensional shapes, flow 
about a modified Shuttle vehicle at angle of attack has been computed at a condition 
for which there exist experimental data as well as numerical results from a different 
computer code. The numerical results were generated using the STEIN code (ref. 7), 
whereas the experimental data (ref. 23) were obtained in the Ames 3.5-Foot Hypersonic 
Wind Tunnel. In figure 14(a), a comparison is shown of numerical and experimental 
surface pressure coefficients on the windward center line plotted against nondimen- 
sional body length for an angle of attack of 25° and a free-stream Mach number of 
10.29. The comparison between all four sets of data is excellent. In fact, the two 
numerical methods give almost the same results. This is interesting since both 
methods have the same resolution in the radial and meridional directions but the 
STEIN code has an axial resolution approximately eight times greater than the HALIS 
code, A further comparison of the experimental and numerical results for this case 
is shown in figure 14(b) where meridional distributions of surface pressure coeffi- 
cients are plotted for three different axial locations on the vehicle. The numerical 
data are only plotted to approximately the tip of the vehicle strake since the geom- 
etry used in the computations has been modified (fig. 5(b)), There is good agreement 
between experimental data and both sets of numerical results for off-axis points. 
Figures 15 to 18 are similar plots of HALIS and wind-tunnel data for angles of attack 
from 30° to 45°. In all figures, there is very good agreement between the experi- 
mental data and the HALIS results. 


Real Gas Computation 

A comparison of perfect-gas (y = 1.4) and real-gas computations for a = 18 

point on a typical Shuttle entry trajectory is shown in figures 19 and 20. A plot of 
center-line pressure coefficients is shown in figure 19. Real-gas effects are most 
prominent in the stagnation region and downstream, where the recompression is weaker 
than that shown for the perfect-gas case. Meridional pressure coefficient distribu- 
tions for both the real- and perfect-gas cases at increasing values of z/L are 
shown in figure 20. 


CONCLUDING REMARKS 

In this paper, a computer code HALIS, designed to compute the three-dimensional 
flow about Shuttle-like configurations at angles of attack greater than 25°, has been 
described in detail. Results from HALIS have been compared where possible with an 
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existing flow-field code; such comparisons have shown excellent agreement. Also, 
HALIS results have been compared with experimental pressure distributions on Shuttle 
models over a wide range of angle of attack. These comparisons have again been 
excellent. It has also been demonstrated that the HALIS code can incorporate equi- 
librium air chemistry in flow-field computations. 


Langley Research Center 

National Aeronautics and Space Administration 
Hampton, VA 23665 
February 10, 1983 
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APPENDIX A 


GENERAL TRANSFORMED EULER EQUATIONS 

In this appendix the components of the Euler equations in both the spherical and 
cylindrical coordinate systems are identified. As well, the development of the 
transformed equations in the respective coordinate systems is presented. 


Spherical System 

In reference 24, the viscous compressible three-dimensional Navier-Stokes equa- 
tions in conservative form are listed for a set of generalized coordinates. From 
these equations with the appropriate coordinates, metric coefficients, and zero vis- 
cosity, the compressible three-dimensional Euler equations in spherical coordinates 
can be written in vector form as 


dW 

5t 


+ 



6G 



Q = 


0 


where w, F, G, H, and Q are defined as 
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(A2b) 


(A2c) 


(A2d) 
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with equations (A3) and the chain rule of partial differentiation, the derivative 
terms in the governing equations can be written in terms of derivatives in the com- 
putational space as 
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Ihe derivatives, found in equations (A4), of the computational coordinates with 
respect to the physical coordinates are determined by using equations (A3) and are 
given as follows; 


an _ 1 

ar r - r, 
s b 


an 

a^ 




a<j> 


(r 




89 




(r 


- 

^b^ae 




at 


at 


8r 


(r - r^) 


b' at 




aa ^ 


8(0 „ 2 

80 It 


(A5) 


27 



APPENDIX A 


By replacing the partial derivatives in equation (A1 ) by their respective representa- 
tion in the transformed space equations (eqs. (A4)), the equations to be solved on 
the computational grid can be written as follows: 


^ 

at a?a^^ a^ a/ 


(A6) 


•The vectors W, F, G, H, and Q are unaffected by the transformation and the 
velocities are still defined in the spherical coordinate system which means that they 
are not necessarily aligned with any grid lines in the computational space. 


Cylindrical System 

From reference 24, the inviscid compressible three-dimensional Euler equations 
in cylindrical coordinates are written as 


where 
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e ^ ^ ^ ~ p{ P/e) . The transformation in the 

cylindrical system is similar to that in the spherical system where now we have 
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Again, the derivatives in the governing equations can be written as follows in terms 
of the derivatives in the computational space: 
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As before, the derivatives in equations (AID) can be determined from equations (A9) 
and are listed as follows: 
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Equation (A7) can be transformed to computational space by replacing the partial 
derivatives with equations (A10) to take on the following form: 


^ + + + + + + + 0 = 0 
5t 5n 5r dh bx] dcj) b<\> dz 6C 


(A12) 


Again^ as in the spherical system, the vectors W, F, H, G, and Q are unaffected 
by the transformation to computational space, and the velocity components are still 
defined in the cylindrical system. 
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EULER EQUATIONS FOR 9 = H 

The Euler equations in the spherical coordinate system, as defined by equa- 
tion (A1), are singular along the ray 0 = te. To obtain a set of equations valid 
along this ray, the following procedure has been used. First, equation (A1 ) is 
differentiated with respect to 9 and then lim 0 ^ u is taken which results in the 
following set of equations: 
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and p = p( p,e) . 


Equations (B1) are valid along a line in space defined by 0 = n. However, this 
line actually represents a spatial plane. In dealing with the flow about 0 = -jt, the 
procedures outlined by Barnwell (ref, 19) are followed. For each value of r along 
the line, the thermodynamic properties are constant for all values of $ as is the 
r component of velocity. However, the 0 and ^ components of velocity will have 
to exhibit some dependence on $ as illustrated in the following sketch: 


- TT 
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In the sketch, we show a 0 surface in the spherical coordinate system with a 
superimposed uniform positive velocity V . Then for any line of constant 0 the 
following conditions apply; 


u , 


V and 

CO 




0 




u . = 0 and u - = V 
0 (}) 00 


(♦= o) 


u 


e 


-V and 
00 




0 


(♦- 



These conditions will always apply along 0 = tu since r 
at this location and lead to the following expressions for 
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The form of equations (B1 ) can be simplified if we choose to solve the equations only 
in the symmetry plane and then distribute the properties about the line 0 = it as a 
function of (|), In the symmetry plane, the velocity component u^ is zero; this 
eliminates the need to integrate equation (Bic) • We can then write the remaining 
equations in the vector form 
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When the derivative definitions (eqs. (A4) ) are substituted into equation (B3)/ the 
following equation, valid in the computational space, is generated: 
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Now r^ and are not functions of $ on 0 = tj. Therefore, 3^1/34) is zero on 

9 = tt; this simplifies equation (B5) to 
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STABILITY ANALYSIS 

M. J, Hamilton 
Sperry Systems Management 
Hampton, Virginia 

Symbols 

A matrix of linearized coefficients of hyperbolic system 

a local speed of sound 

B coefficient matrixes for Euler equations in general coordinates 

Cp specific heat at constant pressure 

c^ specific heat at constant volume 

D weighted superposition of A matrixes used to simplify form of G 

superposition of A matrixes used to simplify form of 
e specific internal energy 

F coefficient functions of hyperbolic system 

F F evaluated at mesh nodes using U 

G amplification matrix for Brailovskaya's scheme 

G^ amplification matrix for MacCormack's scheme 

g^ p abbreviations for quantities appearing frequently in amplification matrix 

for Euler equations in general coordinates 

I identity matrix 

i,j,n,Jl,k indices 

k independent variable in frequency domain 

P similarity transform for Euler equations 

q arbitrary positive integer 

S set of all eigenvalues of D 

t time 

U dependent variable of hyperbolic system 
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approximation to U at mesh nodes from predictor step of Brailovskaya's 
scheme 

U/V,w Cartesian velocities for perfect gas 

V amplitude of Fourier transform component 
Cartesian coordinates 

yj independent variables of hyperbolic system 

a wavelength associated with given wave number and mesh spacing, Ay^ 

Y ratio of specific heats, ^p/^v 

A abbreviation used for terms in amplification matrix 

\ arbitrary eigenvalue of D 

p density of perfect gas 

T right-hand limit of time interval within which entries of G are uniformly 

bounded 

General Analysis 

In this appendix the von Neumann type stability criterion for hyperbolic equa- 
tions in general coordinates is discussed. An explicit linearized stability cri- 
terion is given for application of Brailovskaya's finite-difference method (ref. Cl) 
to the Euler equations. Ihe resultant criterion is also shown to be necessary for 
MacCormack's method (ref. C2) . Finally, the effect of "artificial viscosity" terms 
on the stability is computed. 

Consider the hyperbolic system 
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where U(y^ , . . . t) and Fj^(U) are q-dimensional vectors, and (y-. , . . . ,yj^) repre- 
sents the chosen nonsingular coordinates. Brailovskaya's finite-difference scheme 
for equation (Cl) is a two-step predictor-corrector method with centered differences. 
Applied to equation (Cl), the predictor step is 
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and the corrector step is 
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be treated as constant at each time level. We then apply the usual Von Neumann 
analysis (ref. C3) to the linearized system 
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Let U(k,t) denote the Fourier transform of U(y,t), and let Ay ^ . Set 


^ ^ 


Together, equations (C2) through (C5) give 


n+1 n 

V = GV 


with the amplification matrix 


G = I - D - iD 


where 
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A. sin a. 
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Let the set of eigenvalues of D be denoted by S* From equation (C7), the 
Von Neumann criterion becomes 


\j < 1 for all \ C S 


The scheme of MacCormack does not have the symmetry of the method of Brailovskaya. 
Because of this, the expression for the amplification matrix of MacCormack *s scheme 
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analogous to equation (C7) is not simple. Ihe MacCormack amplification matrix G 
becomes ^ 



2 ^ 2Ai \ «j 

k,j=1 ^ 


1 ) 


+ sin 




sin 




k, j=1 
k<j 


'Vi - 

2 


■[sin 0|^ 


(cos 


a. 

3 


- 1 ) 


1 \ ® 
- sin (cos otj^ - 1 ) J(1 + i) - i At / ^ 

k=1 


sin 


“k 


Ay,, 


(CIO) 


If we assume Oj^ = it, with k = 1, 2, n, the MacCormack problem becomes 

tractable; this is equivalent to assuming alternating signs for the perturbation 
terms at adjacent mesh nodes. In this case, we find 

where 


D 

m 



(C12) 


If S now denotes the eigenvalues of we again find our Von Neumann criterion as 

in equation (C9). 
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Application to Euler Equations 

In Cartesian coordinates, with U = ( p,u,v,w,e)'^, the linearized Euler equations 
for a perfect gas may be written in the form of equation (C4), where 




(Cl 3b) 
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Here p, (u,v,w), and e denote the density, Cartesian velocity, and internal 
energy, respectively. For computational tractability, we diagonalize A 2 * Let a 
denote the local speed of sound, and define P as 


P = 


|~p/e 

0 

0 

0 

1 


0 

1 

0 

0 

0 


0 

0 

0 

1 

0 


PY/a^ 

0 

Y/a 

0 

1 


Then the similarity transform = P^^A^P 


A^ = diag(v, V, V, v+a,v-a) 


For this transformation. 




u 

0 

0 

0 

0 


0 

u 

0 


a2/2y 


-a 2 


a^/2y 


0 

0 

u 

0 

0 


0 
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0 


A 3 = 


w 

0 

0 

0 

0 


0 

w 

0 

0 

0 


0 

0 

w 

a2/2y 

-a^/2y 


0 

0 

Y 

w 

0 


-pY/a^ 

0 

Y/a 

0 

-1 


gives 


0 

"Y 

0 

0 

u 


0 

0 

Y 

0 

w 


(C14) 


(C1 5) 


(C16a) 


(Cl 6b) 
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Now transform to general coordinates. Let (x^,x 2 ,X 3 ) be Cartesian; 

Y ^ ^ nonsingular transformation, in the new coordinates, the linearized 

Euler equation becomes 


5t 


+ 



B. 

D 



0 


(C17) 


wi th 


B. 

D 



+ I 



(C18) 


The results of the section "General Analysis” then hold for equation (Cl 7), with 


D 



B. sin a. 
3 3 



(C19) 


Let 


jA 


formation D = 


denote 

-1 

P DP. 


the Jacobian matrix 5y./5x«, 

y X 

Make the further definitions 


and perform the similarity trans- 


g 


At 



sin a. 
3 




+ va 


3 2 


^ ''“j3 



(C20) 






Aj 




(C21) 
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Combining equations (Cl 5) through (C21) gives 


g 

0 

0 

0 

0 

0 

g 

0 

yPi 

-yPi 

0 

0 

g 

YP3 

-YP3 

0 

a2p^/2y 


^\/ 2 y 

0 

0 

-a2p^/2y 



g - ap 


A straightforward calculation now gives for the eigenvalues of D 


S = 


gf gr g + 



g - 



(C23) 


If the entries of G are uniformly bounded on an interval 0 < At < t, and at 
most one eigenvalue in S has | X] = equation (C9) is sufficient for linearized 

stability (ref, C3) • Since the entries of G are comprised of density, velocity, 
energy, and Jacobian terms, they should be bounded away from coordinate singulari- 
ties; this is assumed to be the case. Then from equation (C23) we arrive at the 
condition 


At 


3 

ua. . + va._ + wa._ + 


j2 33 at 


j=1 


Ay. 


+ a 


3 , V 2 

/a.. 


j,k=1 \ j,k=1 

j<k 


“jl\l ^ “j2\2 °^j3\3 


Ay . Ay, 
3 k 


< 1 (C24) 


For Cartesian coordinates, equation (C24) becomes the well-known 


At 


iHi + IZL + Jwl + al 


1 .1+1 

Ax^ AX3 J 

2 2 2 
Ax^ Ax^ AX3 


< 1 


(C25) 
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In comparing equations (C24) and (C25), we note 


(1) 


Ax 


has been replaced by 


Ay, 




dt 


; thus, stability depends upon the 


contravariant velocity components 


(2) for an orthogonal transformation, the term under the square root in equa- 

.2 


tion (C24) is 


E 


1 


Mesh spacings 


, as in equation (C25) 


For MacCormack's scheme, recall the special case of equation (C1 1 ) • In this 
case, the preceding development holds, with D being replaced by (eq. (Cl 2)). 

One finds equation (C24) again necessary for stability. 

For arbitrary phase, the computational intractability of (eq. (CIO)) makes 

it difficult to establish a general result or worst case. Computational experience 
seems to support the use of equation (C24) for MacCormack*s scheme as well as for the 
more rigorously justified Brailovskaya case. 


Artificial Viscosity Effects 

From equations (CIO) and (Cl 1 ) , it is clear that Brailovskaya's method is not 
dissipative. If sin = 0 for all values of j or for any frequency if D has 
eigenvalues of 0 (in fluid dynamics, this will occur near stagnation or sonic 
points), then G has eigenvalues of 1. This is prevented by adding smoothing or 
damping terms to equation (C6) . A fourth-order term often chosen in fluid mechanics 
is 


n ^ 


-kl 


z 

A=1 






- 4U*? 
D 


+ 6U^ 


JL+^ 


r • • • f , 


/j 


n 


- 4U^ 
D 


/ • • • f _ 






+ 



^A-2' 



(C26) 


where k is the adjustable artificial viscosity coefficient. With this choice, 
equation (CIO) becomes 


G = 


lA - 


iD 


(C27) 
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where 


n 

A = 1 - 2k ^ (cos 2a. 

j=1 ^ 


4 cos a. + 3) 


(C28) 


Equation (C9) now becomes 


( A - + X^l < 1 


(C29) 


This requires 


< 1 


(C30) 


Let X denote the eigenvalues of D as before ^ with X^ At 
tion (C29) gives 


X. Then equa- 



max 


< 


2A - 1 + (5 - 
\| 2 


(C31) 


Our reasoning thus far is general: if equation (C26) is replaced by another 

damping term, equations (C27), (C29), (C30), and (C31) still hold, with equa- 
tion (C28) modified appropriately. Equation (C30) restricts the choice of k. For a 
given k, equation (C31) then restricts our time step. For k = 0, equation (C31) 
reduces to equation (C9) . Equations (C28) and (C30) give the requirement 


0 < k < ^ 


8n 


(C32) 


A restriction 
schemes ( ref . C4 ) . 
appear to be widely 
for the choice (eq. 


of the form of equation (C32) is known for other difference 
Additional restrictions of the form of equation (C31) do not 


recognized. 
(C26)) with 


Figure Cl shows 
n = 2 or 3. 


AtlX 


1 


as a function of k 


max 
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CALCULATION OF PROPERTIES AT SHOCK WAVE 

The shock velocity and other flow properties on the downstream side of the shock 
wave are calculated from the downstream pressure Pg hy using a general three- 
dimensional extension of the axisymmetric method discussed in reference 24* The 
first step in the procedure is to resolve the velocity at the shock wave into an 
orthogonal shock-oriented system with components normal and tangent to the shock wave 
as illustrated in sketch A, 



Sketch A 


The next step is to compute the shock velocity and other flow properties at the shock 
(see sketch A) from the norma 1-shock -wave relations treating the downstream pressure 
Pg and the normal velocity component in the free stream as known quantities. The 
values of p^ used in this step of the procedure are obtained from the finite- 
difference solution of the flow-field equations at the shock. The final step is to 
resolve the velocity at the shock into the components in the coordinate system being 
used to make the flow computations. Since the details of the procedure are dependent 
on the coordinate system employed, the procedure will first be described for the 
spherical coordinate system used in the nose region and then for the cylindrical 
coordinate system used downstream. 
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Spherical Coordinate System 

In spherical coordinates, the shock wave can be described by an equation of the 

form 


F (r,e,$) = r - r (0,^) = 0 
s s 


(Dl) 


where r^CO, ({>) is the radial distance to a discrete set of points which define the 
shock. The outer unit normal to the shock surface is given by the expression 



Thus for spherical coordinates. 


n becomes 
s 


n 

s 


e- 

r 




where 


G 

n 



^r^ sin 0 5$ 



(D2) 


(D3) 


(D4) 


Now a unit vector t^^^ is chosen such that it is tangent to the shock wave and the 
curve formed by the intersection of the shock wave and the plane ^ = Constant. Such 
a unit vector is given by the expression 


t 


si 


ar/d9 
I d5/50| 


(D5) 


where the vector r is defined as follows: 


r = fg(0, $)6j. 


(D6) 
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Thus, 


t 


si 



(D7) 


where 



(D8) 


Now a second unit vector t ^2 tangent to the shock and perpendicular to both Jig 
and tg^ can be obtained from the vector cross product of rig and tg^ 


^s2 


= 


X t 


si 


(D9) 


With equations (D3) and (D7), the following result is obtained for t ^2 



A ^ 

Now these three unit vectors n , t and t ^ defined by equations (D3), (D7), 

s s 1 s2 

and (DIO) can be written in the following short form 


"s = *¥ ®f + Nq6q + 


(D1 1 ) 


^Sl ■ "^ir ®r '^10®e 


(D12) 
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s2 


’’ 2 ? ®r '^29®9 ’^2(J) ®^ 


(D13) 


where 


N- = 1/G 
r n 


(D1 4a) 


% = i7lr)/<=„ 

s 


(D14b) 


1 Sr , 

Nt = -I 3 I^l/G 

lT sin 9 5* 

' s 


n 


(D14c) 
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1r \bQ rt 


(D14d) 


9 = 


(D14e) 




(D14f) 


*’^2r fsin 9 


(D14g) 
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5r 


r sin 9 d<t> 


f(G G. ) 
n t 


(D14h) 


T - = 
2 (|) 


r + 
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(D14i) 
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Note that for the special case of Q = n and ^ = ±n/2r N-f T T / and G 
reduce to the following limiting forms: (J) 2 r 20 ^ 


N-=T- = T ^=0 
4 ) 2 r 20 






1 


5r ^ 


(0 = 71 and = ±%/2) 


1 + — 


- 90 

r 

s 


(D15) 


The free-stream velocity in spherical coordinates is 


V = V- e- + V. e. + V- e- 
® r,«r 0,®0 


(D16) 


where 


V- = |v I (sin a sin 0 sin * + cos a cos 0) 

^ 00 ' CD ' 

== I (sin a cos 0 sin A - cos a sin 0 ) 

9 , 00 ' oo' ^ 


V- = Iv 

^-00 I 00 I 


sin a cos 


(D17a) 

(D17b) 

(D17c) 


The free-stream velocity component normal to the shock wave V is 

n, “ 


V = V • n = N- V- + N V + N- V- 
n,-» ® s r r,oo ^ "e''e,oo ^ 


(D18) 


and the free-stream velocity components tangent to the shock wave V and V 

are t 2 ,» 


V = V • t 

t 1 , oo 00 si 


T-v- +T V +T-V- 
1r r,cn ^ 10 '' 0 ,oo 1 (J) ^4), 00 


(D19) 


and 


=V *t^ = T^~ V- + T V + T - V- 
t2,«> CO s2 2r r,oD 20 0,® 24) 4>,c 


(D20) 
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Although, the tangential components of velocity and v^2 unchanged across 

the shock wave, the other properties illustrated in sketch B change and must be 



Sketch B 


computed from the following norma 1-shock -wave relations for a moving shock: 

p/v -v\=pfv -v\ 

n, «> sy n,s sj 


p + pfv -vV=P + p/v *-vV 

® n,® sj s ^s I n,s sJ 


(D22) 


h + 

00 


(v - V (v - V 

\n,® s/ , .\n,s s/ 

2 = ^ 


(D23) 
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or, since = 1 / equations (D21) and (D22) become 



Now, rearranging equation (D25) results in the following equation: 


(D24 

(D25 


P 


s 


(V 

\ n,® 



fv 

\ n,s 


(v 

\ n,® 



(D26 


Combining this equation with equation (D24) and solving for (v - V gives the 
following equation: ® 


2 Ps - P® 


fw - V \ = - 

V n,® s) 1-1 


/Pc 


Similarly, equation (D23) becomes 


(D27 



Equations (D27) and (D28) can be combined with an equation of state in the form 
p = pidi 

to yield the following result: 

i_ _ 2h^ + (p^ - P^) 

PS ' (2Pg/K) - (Pg - P„) 


(D28 


(D29 


(D30: 


where for an ideal gas 
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and for a real gas, < is evaluated from p, p, and h at the previous time step 


K = 



(D32) 


Equation (D32) is a good approximation of the value of k at the early time steps 
and it becomes more and more "exact" as the solution approaches convergence. This 
process eliminates the need to iterate the shock-wave relations by using the "real- 
gas subroutines" and, thus, speeds up the solution considerably. With equa- 
tion (D30), equation (D27) can be solved to obtain the following equation for the 
shock velocity Vg-. 



(D33) 


Then use this result in equation (D24) to obtain ^ as 

ri , s 


V = —(V - V ) + V (D34) 

n,s n,® s s 


The enthalpy h^ can then be computed from equation (D28) and e^ from the relation 


e = h 
s s 



(D35) 


The free-stream velocity components in the spherical coordinates V- , V- _ f 

r , ® u 

and V- can be transformed to shock -oriented velocity components v_ _/ V.- r 

(|),oo n,oo t1 ,® 


and V ^ by using the transformation matrix A 
t2, <» 


1 1 

®12 

®13 


^1r 

^10 

^1^ 


21 

®22 

®23 

= 

2r 

'^29 


(D36) 
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as follows: 


V- 



r, ® 


t1 ,«> 

V. 



V ^ 

8 

CD 


t2,<» 

V-. 


V 

L 


n, «> 


(D37) 


This same transformation matrix can be applied to conditions behind the shock wave; 
thus, 


u- 



r,s 


t1 ,s 



V ^ 

9,s 


t2,s 



V 

- 


n,s 


where 


(D38) 


V = V 

t1,s t1,< 


V = V 

t2,s t2,< 


(D39a) 

(D39b) 


TSie spherical velocity components downstream of the shock wave can be determined from 
equation (D38) as follows: 


u- 




r,s 


t1 ,s 



= A-^ 

V 

(D40) 

0,s 


t2,s 


u- 


V 


(j),S 


n, s 
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where is the inverse of the transformation matrix A* Now can be written 

in the general form 


\ I Adj A 1^1 


1 1 

Si 

‘12 

S 2 

^13 

S 3 


31 


32 


'33 


(D41 ) 


where |a| is the determinant of the matrix A and A^j is the cofactor of the 
element j • ^ spherical coordinate system we can write 


‘I = - ^ 2 * “ 9 ) - ’’i9(’^2; 


(D42) 


and 


^11 '^2e”i ” ’’^2^ ^0 


A = -T - N- + T^- N- 
12 2r (|) 2(1) r 
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Now from equations (D40) and (D41), the spherical velocity components downstream of 
the shock are obtained as follows: 


'^0,s ~ |a|(^12\i,s * ^22^t2,s * ^32\,s) 
'^4),s ^ ,s ^23^2,8 S3^n,s) 


(D44) 


(D45) 


(D46) 


Cylindrical Coordinate System 

In cylindrical coordinates, the shock wave can be described by an equation of 
the form 


F (r,(J), z) =r -r ((j>, z) 
s s 


(D47) 


In a procedure similar to that used with the spherical coordinate system, three unit 
vectors related to the shock wave are defined. 

The first, ng , is the outer unit vector to the shock wave and is defined by 
ng = + N0eQ + N^e^ (D48) 

where 
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'/°n 


(D49a) 


N 





(D49b) 


N = 
z 



(D49c) 


60 


APPENDIX D 


and 


G 

n 





(D50) 


The second, unit vector tangent to the shock wave and the curve 

formed by the intersection of the shock and the plane 4> = Constant and is defined 
by 


A 

"^sl 

= TirSr + Ti4,e^ + T-^z&z 

(D51) 


/ar \ / 


1r 
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(D52a) 

T 

= 0 

(D52b) 

Iz 

= 1/G^ 

(D52c) 


and 



(D53) 


The third, £ 32 ^ is a unit vector perpendicular to both £ and toi and is 
defined by ' 

^s2 “ "^ 2 r®r '^2<|)®4) ■*" *^ 22^2 (D54) 

where 
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(D55a) 
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^ 24 >= ^ 


f(G G ) 
n t 


^2z = 


The free-stream velocity in cylindrical coordinates is 


V=v e+V e+V e 
00 r , CO r <(),<» (j) z /oo z 


= V sin a sin 


= V sin a cos (J) 


= V cos a 


The free-stream velocity component normal to the shock wave is given by 


V = V • n 

n, oo 00 s 


N V + N V + N V 
r r,oo 0^00 z z,oo 


and the velocity components tangent to the shock wave are given by 


= V • t , = T V + T, V, + T V 
t1,oo 00 si r,oo 1 <|) 4),oo 1z Z ,< 
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t2,oo 00 s2 2r r,oo 2<t) 2 z z ,<x> 

Now the tangential components of velocity are unaffected by the shock wave; thus, 


V = V 

t1,s t1,< 
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The velocity ^ and other properties on the downstream side of the shock wave 
(sketch B) are c6mputed by using equations (D23) through (D35) in exactly the same 
way as for the spherical coordinate system described previously* 

Similarly, the cylindrical velocity components on the downstream side of the 
shock wave are computed from the equations 


1 , 

r,s |A| ^ 11 t1,s 21 t2,s 
1 , 

<t>,s |A|\ 12 t1,s 22 t2,s 
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z,s |a| \ 13 t1,s 23 t2,s 

where 

|a| = T, N - N \ + T 

' ' 1r^ 24 ) z 2z 4>j 1z 
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(a) Complete geometry. 


(b) Modified geometry. 


Figure 5,- QUICK geometry models. 
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Profile view 


Cross section at z = 50 in 



Cross section at z = 350 in. Cross section at z = 650 in 

Figure 8.- Initial shock shapes, = 10*3; a = 25° 
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Figure 13.- Surface axial Mach number distribution 








(a) Windward center-line pressure coefficients. 
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(b) Meridional pressure coefficients. 

Figure 14.- Comparison of experimental and calculated pressure distributions at 

M_ = 10.29 and a = 25». 
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(a) Windward center-line pressure coefficients. 

Figure 15.- Comparison of experimental and calculated pressure distributions at 

= 10.29 and a = 30®. 
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(a) Windward center-line pressure coefficients* 
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(b) Meridional pressure coefficients. 


Figure 16.- Comparison of experimental and calculated pressure distributions at 

= 10.29 and a = 35® . 
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